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ABSTRACT 

We introduce and analyze several aspects of a new model for cell differentiation. 
It assumes that differentiation of progenitor cells is a continuous process. From the 
mathematical point of view, it is based on partial differential equations of transport 
type. Specifically, it consists of a structured population equation with a nonlinear 
feedback loop. This models the signaling process due to cytokines, which regulate 
the differentiation and proliferation process. We compare the continuous model to 
its discrete counterpart, a multi-compartmental model of a discrete collection of cell 
subpopulations recently proposed by Marciniak-Czochra et al. [T7j to investigate the 
dynamics of the hematopoietic system. We obtain uniform bounds for the solutions, 
characterize steady state solutions, and analyze their linearized stability. We show 
how persistence or extinction might occur according to values of parameters that 
characterize the stem cells self-renewal. We also perform numerical simulations and 
discuss the qualitative behavior of the continuous model vis a vis the discrete one. 
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Introduction 

Cell differentiation is a process by which dividing cells become specialized and equipped 
to perform specific functions such as nerve cell communication or muscle contraction. Dif- 
ferentiation occurs many times during the development of a multicellular organism as the 
organism changes from a single zygote to a complex system with cells of diff'erent types. 
Differentiation is also a common process in adult tissues. During tissue repair and dur- 
ing normal cell turnover a steady supply of somatic cells is ensured by proliferation of 
corresponding adult stem cells, which retain the capability for self-renewal. Also various 
cancers are likely to originate from a population of cancer stem cells that have properties 
comparable to those of stem cells |3] . 

Stem cell state and fate depends on the environment, which ensures that the critical 
stem cell character and activity in homeostasis is conserved, and that repair and develop- 
ment are accomplished [I9]. Cell differentiation and the maintenance of self-renewal are 
intrinsically complex processes requiring the coordinated dynamic expression of hundreds 
of genes and proteins in response to external signaling. During differentiation, certain 
genes become activated and other genes inactivated in an intricately regulated fashion. 
As a result, differentiated cells develop specific structures and performs specific functions. 
There exists evidence that disorder in self-renewal behavior may lead to neoplasia [31 H]. 
For example, it has been shown that acute myeloid leukemia originates from a hierarchy of 
cells that differ with respect to self-renewal capacities [T3|[23]. Although much progress has 
been made in identifying the specific factors and genes responsible for stem cells decisions 
|2U] . the mechanisms involved in these processes remain largely unknown. 

While different genetic and epigenetic processes are involved in formation and mainte- 
nance of different tissues, the dynamics of population depends on the relative importance 
of symmetric and asymmetric cell divisions, cell differentiation and death. The same genes 
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and proteins are observed to be essential for regulation of different tissues [21]. This unity 
and conservation of basic processes implies that their mathematical models can apply 
across the spectrum of normal and pathological (cancer stem cells) development. 

One established method of modeling such systems is to use a discrete collection of ordi- 
nary differential equations describing dynamics of cells at different maturation stages and 
transition between the stages. These so called multi-compartmental models are based on 
the assumption that in each lineage of cell precursors there exists a discrete chain of matu- 
ration stages, which are sequentially traversed, e.g., [151 [27]. However, it is also becoming 
progressively clear that the differentiated precursors form such sequence only under home- 
ostatic (steady-state) conditions. Committed cells generally form a continuous sequence, 
which may involve incremental stages, part of which may be reversible. As an example, 
cell differentiation without cell divisions is observed during neurogenesis. Moreover, in 
some tissues such as the mammary gland, different stages of differentiation are not well 
identified [5]. 

These observations invoke not only the fundamental biological question of whether 
cell differentiation is a discrete or a continuous process and what is the measure of cell 
differentiation, but also how to choose an appropriate modeling approach. Is the pace of 
maturation (commitment) dictated by successive divisions, or is maturation a continuous 
process decoupled from proliferation? In leukemias, it seems to be decoupled. The classical 
view in normal hematopoiesis seems to be opposite. 

To address these questions and to investigate the impact of possible continuous trans- 
formations on the differentiation process, we introduce a new model based on partial differ- 
ential equations of transport type and compare this model to its discrete counterpart. The 
point of departure is a multi-compartmental model of a discrete collection of cell subpop- 
ulations, which was recently proposed in [T7] to investigate dynamics of the hematopoietic 
system with cell proliferation and differentiation regulated by a nonlinear feedback loop. 
Furthermore, since self-renewal is an important parameter in our models, the proposed 
models seem to be a right departure point to investigate cancer development, for example 
in leukemias [23] . 

In the present paper we extend the discrete model to a structured population model 
accounting for a continuous process of differentiation of progenitor cells. Models of the 
latter type have been already applied to the description of some aspects of hematopoiesis 
[HI |5l [H |2l [11] . These models are based on the assumption that differentiation of progenitor 
cells is a continuous process, which progresses with a constant velocity. Mathematical 
description involves so called age-structured population equations. The model presented 
here is novel due to the nonlinearities in the coupling of the model equations, in particular 
the nonlinear coupling in the maturity rate function. 

The paper is organized as follows. In Section [T] we formulate the new model. In Section 
ll.2l the link between our model and the discrete model of [17] is accomplished. Sections [2l 
[3] and [H are devoted to the analysis of the model. In Section [2] the existence and uniform 
boundedness of the solutions are shown. In Section [31 it is shown that depending on the 
value of a parameter characterizing stem cells self-renewal model solutions tend to zero 
or they stay separated from zero. Section [4] provides the structure of steady states and 
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conditions for existence of a positive stationary solution, while Section 14.21 is devoted to 
a linearized problem around the positive steady state when this state exists to investigate 
its stability. Using the characteristic equation we study some special cases, for which we 
show stability or instability of the positive stationary solution. In Section [5l a numerical 
approach and some results on stability and instability are presented. We conclude in 
Section [6] with some final comments and suggestions for further investigation. 



1 Model of Cell Differentiation 
1.1 Continuous Model 

In the following we assume that the dynamics of differentiated precursors can be ap- 
proximated by a continuous maturation model. Under this assumption we extend the 
multi-compartmental system from [17]. Let w{t) denote the number of stem cells, v{t) the 
number of mature cells and u{x, t) the distribution density of progenitor cells structured 
with respect to the maturity level x, so that J^^ u{x, t)dx is equal to the number of progen- 
itors with maturity between xi and X2- This includes maturity stages between stem cells 
and differentiated cells. Thus, M(0,t) describes a population of stem cells and t), for 
X > 0, corresponds to progenitor cells. We assume that x = x* denotes the last maturity 
level of immature cells, and therefore, u{x*,t) describes the concentration of cells which 
differentiate into mature cells. 
The model takes the form 

[2a^(s) - l]pw{s)w{t) - dyow(t), (1) 

p{x, s)u{x,t) — d{x)u{x,t), (2) 
2[l-a^,{s)]p^{s)wit), t>0, (3) 

g{x*,s)u{x*,t) - fivit), (4) 

together with initial data 

w{0) = wo>0, n(0, x) = uo{x) > 0, v{0) = vo>0. 

Integrating formally Equation ([2]) and adding it to Equations ([1]) and (jl]) yields the follow- 
ing cell number balance equation 

d r* r* 

— \w{t)+ I u{x,t)dx + v{t)^ = ipyjis) — dw)w{t) + I (p(x, s) — d(x))u(x,t)dx — fiv(t). 
dt' Jo io 

(5) 

System describes the following scenario: After division a stem cell gives rise to two 

progeny cells. Cell divisions can be symmetric or asymmetric. We assume that on the 
average the fraction of progeny cells remains at the same stage of differentiation as 
the parent cell, while the 1 — a^j fraction of the progeny cells differentiates, i.e. transfers 



d_ 

dt 



w{t) 



dtu{x, t) + d^[g{x, s)u{x, t)] 
g{0,s)u{0,t) 

d . , 
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to the higher differentiation stage. This covers the symmetric and asymmetric scenarios. 
Parameters and dw denote the prohferation rate of stem cells and their death rate, 
respectively. Progenitor cells differentiate at the rate g, which depends on their maturity 
stage and is also regulated by the feedback from mature cells given by a signaling factor s. 
Parameters p{x) and d{x) denote the proliferation and death rates of precursor cells and 
depend on the level of cell maturation. Mature cells do not divide and die at the rate /i. The 
whole process is regulated by a single feedback mechanism based on the assumption that 
there exist signaling molecules (cytokines) which regulate the differentiation or proliferation 
process. The intensity of the signal depends on the level of mature cells, and is modeled 
using the dependence 

s = s[vit)] ^ 



which can be justified using a quasi-steady state approximation of the plausible dynam- 
ics of the cytokine molecules, see [17]. This expression reflects the heuristic assumption 
that signal intensity achieves its maximum under absence of mature cells and decreases 
asymptotically to zero if level of mature cells increases. 

The concentration of signaling molecules s{v) influences the length of the cell cycle 
(proliferation rate p) and/or the fraction of stem cells self-renewal (a^) as well as the rate 
of cell maturation (g). 

In this model, differentiation of stem cells takes place during mitosis. The differentiation 
of progenitor cells occurs independently of proliferation. In other words, cells undergo 
continuous transformations between divisions. We call this process maturation. In the 
terms of the model this means that in an infinitesimal time interval (t, t + dt), the following 
events occur to a cell of maturity x, 

1. either the cell matures to level x + dx, which happens with probability g[x,v(t)]dt, 

2. or the cell divides into 2 daughters, which happens with probability p{x)dt, with 
other events occurring with probabilities of the order o{dt). 

If we stick to the discrete model proposed in [17], we obtain the following relations 
linking proliferation and maturation 

g{x,v) = 2[l-j^]p{x), 

a«, = a(0), Pw=p{0), < a^„ = a(0) < 1, 

which does not necessarily mean that differentiation can only occur by division; see the 
discussion at the end of Section II. 2[ 

This leads to a simplification of the boundary condition (|3]) that becomes 



u{0,t) = w{t). 
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1.2 Discrete versus Continuous Models 

General Setting 

In this section we consider the relationship between the structured population model ([T])- 
(jl]) and the multicompartmental model introduced in [17]. Following reference [18], the 
multicompartmental model can be formulated in a general way, 

Pi{s)ui - gi{ui,s) - diui, (7) 

Pi{s)ui + gi^i{ui_i,s) - gi{ui,s) - diUi, for i = 1, n - 1 (8) 

5f„_i(u„„i, s) - d„,u„, (9) 

where gi{ui, s) denotes a flux of cells from the subpopulation i differentiating to the sub- 
population i + 1. The terms Pi{s)ui and diUi describe cell fluxes due to proliferation and 
death, respectively. In the general case proliferation or differentiation may depend on 
signal intensity. 

In reference [17] , differentiation was linked to proliferation and the following expressions 
were proposed: 

gi{ui, s) = 2(1 - ^ ^ )PiUi. (10) 

This is similar to what was defined for stem cells with p^ and a^; here represents 
the fraction of cells remaining at the same stage of differentiation i as the parent cell while 
the 1 — fraction of cells differentiates to the higher stage i + 1. Formulation ([7])- ([2]) 

describes the differentiation process independently of cell proliferation in the sense that 
cells either multiply at stage i or differentiate from compartment i to i + 1 and so forth. 
Assuming that cell differentiation occurs at a properly-chosen time scale compared to the 
time scale of the cell division process, we show in the next paragraph how to obtain, after 
a suitable renormalization, the structured population model of the next paragraph. 

Continuous Limit 

Let us write System ([7])-(|n]) in a dimensionless way. We define V, V, Qi, Q, Ui, hi and W„ 
as characteristic values for the quantities pi, di, gi, gi for i > 2, ui, Ui for 2 < i < n — 1, 
and Un, respectively. Then, we define dimensionless quantities by = ^ etc. We make 
the following hypothesis 

gi{ui, s) = gi{s)ui. 



d 



dt 
d 

—u. 
dt 

d 
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Then, system ([7])-([9]) becomes 



= Vpi{s)ui-gigiUi-VdiUi, (11) 

dt 

—U2 = Vp2{s)u2 + Ql-^giUi - Gg2U2 - Vd2U2, (12) 

dt U 

—Ui = Vpi{s)ui + Q{gi^iUi^i - giUi) - VdiUi, for i = 3, ...,n - 1, (13) 
dt 

-^Un = —gn-lUn-1 ~ VdnUn- (14) 

Letting the number of compartments tend to infinity, we pass from the discrete model 
to the continuous model by associating to the Wj's a function, constant on intervals of 
type {ei,e{i + 1)), with 5 — )■ 0, i — )■ oo and the product ei remaining positive and finite, 
say n = n^, with en — ?■ x* G M.*^ = (0, +00). Compartment dependent constants tend to 
continuous functions, sums over the index i are interpreted as Riemann sums tending to 
integrals while finite differences give rise to derivatives. A precise discussion of the limiting 
process is outside the scope of this presentation. We refer, for instance, to [3 [10] for recent 
examples of how to obtain such limits based on moments estimates. 

In order to interpret the terms Q {gi-iUi^i — giUi) in Equation 0131] and Gijj-giUi — Qg2U2 
in Equation f[T^ . as a finite differences tending to a derivative, we take 

^1 ^ Wi 1 

e U e 

Assuming that Vpi{s)ui and VdiUi tend toward limits p{x, s)u{x,t) and d{x)u{x,t) leads 
to 

V = l, V = l. 

In order to obtain that jj^gn~iUn~i converges to the limit g{x*)u{t, x*) in Equation (111]), 
we require 

^ _gu _iu _u 

~ TT ~ ~tT " ~ ~' 

This means that the order of magnitude of the number of mature cells is much larger than 
the one of the maturing cells and stem cells. This can be interpreted by the fact that Ui 
tends to u{t,x) a density of cells per unit of maturity, whereas Ui and m„ are numbers of 
cells. This scaling follows also from mass balance considerations: System (I7])-(|n]) leads to 
the following mass balance 

d d n—i 71—1 

^ X] ^ P^^^ X^Pi^ii - diUi - ^ diUi - dnUn, (15) 

i=2 1=2 i=2 

meaning that the exchange among compartments at maturation rate gi does not influence 
the total growth of the population. In this equation, we keep the specific values Ui and m„ 
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and interpret the sum for 2 < i < n — 1 as an integral. We obtain 

i=2 i=2 

which leads to the following choice 

U _U _ 

With this choice and the previous relations, we are led to choose Qi = 1, which allows a 
limit for Equation (ITTl) . We note however that this implies a different order of magnitude 
for Qi and for Q; the interpretation could be that we have divided the previous discrete 
compartments into smaller ones, of size e, where division does not occur but where matu- 
ration occurs. In this framework, Qi is not homogeneous to Q but rather to the integral of 
Q over a small compartment of size e. 
Under these assumptions, let us set 

Xt{x) = Xlie,{i+l)e){x), 

with xa being the indicator function of a set A. We introduce the piecewise constant 
function 

u%x,t) := ^Mi(t)x-(a;). 
By the same token, we associate the following functions to the coefficients 

rig rij rij 

d%x) := ^dixtix), p%x,s) := (x), g'{x,u{x),s) := ^ gi{ui, s)xt{x). 

i=l i=l i=l 

We make the following continuity assumptions on the dimensionless system: 

3K > s.t. \gi\ + \di\ + \pi\ < K, \gi+i - gi\ + \di+i - di\ + \pi+i - Pi\ < f 

(16) 

Pi,gi are uniformly continuous with respect to the variable s. 

We define the piecewise constant functions g^, d'^ and p'^ on the respective basis of the 
discrete coefficients gi, di and Pi, similarly as u'^ was defined for Uj. Assumption f lT6|) leads 
to their convergence (up to subsequences) to continuous functions g, d and p of both 
variables x and s (see Lemma 1 of [TU] for instance). We can prove (based e.g. on [TJ [TU] ) 
the following result. 

Proposition 1.1 Suppose that is a solution of Equation f|T3|) verifying {ui{t = 0)) G l^. 
Under Assumption f lT6|) . for allT > 0, there exists a subsequence of [uf) converging towards 
a limit u G C{0,T; Ai^{[0,x*]) — weak—*) solution of Equation (j2]), ul to a limit Ui G C(0,T) 
solution of Equation ([T]) where s = s{v) with v a limit of a subsequence ofu^^. The boundary 
condition (E]) is satisfied in a distributional sense and Equation ^ is satisfied, what is 
equivalent to a weak formulation of the boundary condition (|4]) . 



A structured population model of cell differentiation 



9 



We note that we have used the fact that maturation and prohferation are decorrelated. 
If this were not the case, it would be impossible to make a 1/e factor appear in Q, since 
V = 1. In such case in the limit equation the transport appears as a first order corrective 
term and Equation (Ej) is replaced by 

dtu{x,t) + edx[g{x, s)u{x,t)] = p{x, s)u{x,t) — d{x)u{x,t). (17) 

Figure [T] in Section 15.21 depicts related numerical simulations. 

2 Uniform Bounds for the Continuous Model 

In the remainder of this work we will consider a special version of the above model assuming 
time independent proliferation rates p{x), and zero death rates of undifferentiated cells 
= and d{x) = 0. Indeed, neglecting death rates of immature cells does not change the 
analysis. Concerning the feedback loops it was shown in [T7] for the discrete model that 
the feedback on the stem cells self-renewal fraction and on the maturation speed g is much 
more important for the efficiency of the process than the feedback on the proliferation 
rate p{x). Therefore, in the reminder of this work we focus on the model with regulated 
self-renewal and maturation. We also introduce simpler notation which makes it easier for 
analysis. This yields the following system of differential equations for t > 0, a; > 0. 



-w{t) = a{vit))w{t), (18) 

dtu{x,t) + d^[g{x,v(t))u{x,t)] = p{x)u{x,t), (19) 

u{0,t) = w{t), (20) 

j^vit) = gix*,v{t))uix*,t)-fxvit), (21) 

together with initial data 

w(0) =wo>0, u{0, x) = uoix) > 0, v{0) = vo>0. (22) 
We obtain the cell number balance law 

— [u; + f u{x,t)dx + v] = [a{v) + g{0,v)]w + f p{x)u{x,t)dx — nv{t), (23) 



corresponding to the fact that the total population can only change by proliferation or 
death. Indeed, one can interpret a + g{0,v) as the stem cells proliferation rate, see above 
sections and Equation ([5]). 

In the sequel we will study model f|T8|) - fl2T]) under the following assumptions 

g., g.. G L°°([0,a;*] X R+), a{x) E C([0,cx))), pix) E C\[0,x*]), (24) 

a{v) E [aoc^o], « is decreasing , a(+oo) := Ooo < 0, (25) 
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< 5- < g{x,v) <g+<oo, V(x,t;) G [0,x*] x M+, (26) 
First we show that the model solutions are uniformly bounded 

Theorem 2.1 Under assumptions (!24|) - (!26l) and that Uq{x) G C^([0,x*]), the solution to 
System (fT8|) -( 122|) uniformly hounded. More precisely, all the components w{t), u{x,t), 
v{t) are uniformly bounded. 

The remainder of the section is devoted to the proof of this result, which uses some technical 
lemmas. We first prove the following estimate 

Lemma 2.2 Under the assumptions of Theorem \2.1\ the function z{x,t) = dx{hiu) is 
uniformly bounded on [0,x*] x ]R+. 

Proof. The equation for z reads 

dtz + dx{gz) = -g^x + p^, 

^(0,^) =-^^-^eL-(0,+oo). 

Indeed, we have z{0,t) = ^^ttty' ^^^^ we can compute 



z(0,t) 



u{0,t) 

dtujO, t) - p(0)m(0, t) + g,{0, v)u{0, t) 

giO,v)uiO,t) 
a{v)-p{0) gx{0,v) 



g{0,v) g{0,v)- 



And we conclude that z{0,t) is uniformly bounded by assumptions ( l24l) - (!26l) . 
Next, we rewrite the equation for z as 

dtZ + gd.^z = -g^z + Q{x,v), (28) 

where Q = —gxx + is a bounded function of v and x. 

In the following, we show that the solution to ( l28l) satisfies the estimate 

\\z{xM\l^ < M := ( sup|^(0,t)| +sup|2(a;,0)| +x* ||-||l- | e^*"tll^°°. (29) 
\ t X 9 J 



Indeed, since 5^ > > 0, we can rewrite Equation ( 1281) as 

g g g 

and apply the method of characteristics by defining as usual (except that x plays the role 
of time), 

diT 1 

— (x, t) = -(x, v{T{x, t))), T{x = 0,t)= t, 
ax g 
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-^(x,t') = —(x,v(T-\x,t'))), T~\x = 0,t')=t'. 
dx g 

We look for solutions of the form Z{x,t) = z{x,T{x,t)), which satisfy the following equa- 
tion 

d.Z + ^{x,T{x,t))Z = d^{zJ '''' '''''')e '''"=^{x,v{T{x,t))). 
9 9 

Integrating the above equation yields 

J 9 



Z{x,t) = Z{0,t)e +j ^(e,^(T(e,t)))e « d^, 



-/^(C,T(c,i))dC }o -/^{C,T(c,t)K 
^(x,T(x,t)) = ziO,t)e ^ + ^(e,^(r(e,t)))e « d^. 

J 9 



Defining t = T{x, t), or yet t = T~^{x,i), yields > for i > and we obtain 



z{x,t) = z{0,T-\x,t))e ' ' ' + ^(e,^(r(e,T-i(x,t))))e ^ d^ 

J 9 





Therefore, for t > — — it holds 



\z\\L^<{\z{S).-)\+x*\\^\\L^)e^ I'T 

9 



From Lemma [2.21 we deduce several useful estimates 



Lemma 2.3 There exist positive constants Mi, M2, M3 such that the solutions to system 
(fT8|)-(|22|) satisfy 



(i) w{t) < Miu{x,t), 

(ii) w{t) < M2v{t), 
(Hi) u{x,t) < M-iw{t). 
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Proof (i) Boundedness of —z = — ^ Inu results in the following inequality 

In - < In — + Mx, 

u w 

which in turn yields assertion (i) with Mi = e^'^* . 
{ii) To bound w by f , we calculate 

at V V \ V 

Since a{v) < a{0), g{x*,v{t)) > and u{x*,t) > w(t)/Mi we obtain 

d w w f g_ w\ 

dt V V \ Ml V J 

This yields the estimate 

Mt)<.(t)maxf^,M,^^i^) :=MMt). 
Vf(0) 9- J 

and the assertion {ii) is proved. 

[in) The proof follows as in (i), departing from \n.u{x,t) < ln'u;(a;,t) + Mx. □ 

As a consequence of Lemma 12.31 we derive 

Corollary 2.4 Under the assumptions of Theorem \2.1\ the components w{t), u{x,t) and 
v{t) of the solutions to System fITS]) - fl^2]) are uniformly bounded. 

Proof. Applying Lemma [2.31 (ii) to equation f|T8l) . we obtain 

dw ( w , 

— ;— < tt W. 



dt ~ \M: 



'2 



This yields boundedness of w by Assumption ( 1251) . 

Boundedness oiw yields also boundedness of u using Lemma [2.31 iiii). Finally, bound- 
edness of V results from Equation f[2Tl) due to boundedness of u{x*,t) because g < g+. □ 

The proof of Theorem 12.11 is now complete. □ 



We also state another result, in the spirit of Lemma [2. 3[ that is used later on 
Lemma 2.5 There exists a constant M4 > and < 7 < 1 such that v{t) < M^w^'lt). 
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Proof. We calculate 



d V 
dt w'l 



<Msg{x*,v{t))w 



r(/^ + 7«(^))- 



We choose 7 > small enough such that /x + 7000 '■= Hi > and 7 < 1. Since w is 
uniformly bounded, we find 



Finally, we conclude this section with a consequence of Theorem 12.11 

Corollary 2.6 Under the Assumptions (124|) -(l26l) andthatuo{x) G C^{[0,x*]), System (fT8|) - 
(!22|) has a unique global solution. Furthermore, such solution is uniformly bounded. 

Proof. Local in time existence of the unique solution follows from the Cauchy-Lipschitz 
theorem. Theorem 12.11 provides uniform boundedness of solutions and hence the global 
existence. □ 



3 Extinction and Persistence 

In this section we provide conditions for extinction and persistence of positive solutions. 

First, we consider a case when a(0) < 0. In this case there exists only a trivial steady 
state of the model and 

Theorem 3.1 Assume (|24D - (|26|) . // a(0) < 0, then all solutions of system (fTSD - (122|) 

converge to zero at an exponential rate. 

Proof First of all, notice that, since a{v) < a(0) < 0, it is obvious from equation f fTSj) 
that w converges to exponentially. 

For the other components, we consider a functional 'jw(t) + e~^^u{x,t)dx + e~'^^*v, 
with positive constants 7 and /3 to be determined. We compute its time derivative. 




which yields boundedness of 



a 




V 
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Since a{v) < a(0) < we may choose 7 such that 7a(0)+sup„ ^'(0, v) < —T < 0. Moreover, 
choosing /3 such that p{x) — f3g{x,v) < — F imphes that 



d 
dt 



{-fw{t) + / e-^^M(x, t)dx + e-^'^'v) < -Tw{t) -T e-^^M(x, t)dx - ^le-^'^'v. 
Jo Jo 

We conclude that the solutions converge to zero at an exponential rate for t — > 00. 



Secondly, if it is the case that a(0) > 0, then we conclude that the solutions to the 
system cannot become extinct. 

Theorem 3.2 Assume (l2l) -fl26l) . w{0) > and uo{x) G C\[0,x*]). If a{0) > 0, the 
solution u, V, w of system (IT8l)-( !2T]) with positive initial conditions remain hounded away 
from zero. 

Proof Applying Lemma [2.51 to equation ( lT8l) . we obtain 

dt 

and the assumption a(0) > allows us to conclude. Then, the estimates of Lemma [2.31 
conclude for u and v. □ 



4 Stationary Solutions and Their Stability 
4.1 Stationary Solutions 

As usual in dynamical systems, a natural question concerns the existence of steady states 
(stationary solutions). We shall now investigate this issue. 

In our case, the steady states are given by the solutions {w, u, v) to the system 

a{v)w = 0, (30) 
-^[g{x)u{x)] = p{x)u{x), (31) 

m(0) = w, (32) 
g{x*)u{x*) - fiv = 0, (33) 

where g{x) := g{x,v). 

System (ITS]) - ([^T]) always admits the trivial steady state w = 0, u = 0, v = 0, which 
we do not consider. Depending upon the value a(0) it may also have exactly one positive 
steady state [w, u, v), as we state it in the 

Lemma 4.1 Under the Assumptions 0241] - (12^ . the System 0181] - 02 ip has a strictly positive 
steady state if and only if a{0) > 0. Furthermore, the steady state is unique. 
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This condition is in agreement with biological observations concerning self-renewal of 
stem cell subpopulation [12] and an analogous condition for the compartmental model was 
discussed in [T8l. 



Proof Since we discard the trivial steady state, from equation ( I30l) we obtain the con- 
dition a{v) = 0. As we know that a decreases and tends to Ooo < at infinity, there exists 
a unique solution v to 

a{v) = 0, (34) 
if and only if the condition a(0) > holds. Thus, we may compute 

u{x ) = — — - . 

g{x*} 

Solving differential equation ( 13 ip with this boundary condition a,t x = x* yields, 

= ^uix*)exp { - r ^ciej. (35) 

We finally identify iv using the boundary condition at x = 0, iv = u{0), which leads to 

w = ——u{x ) exp \ - -TT^dn = —— exp \ - ^ttt^Q- 36 
^(0) L Jo giO J ^(0) L Jq g{^) J 

This gives explicit values of the model and completes the proof. p 
For a given explicitly by (E]), we may compute 



u[x 



2am — 1 
k ' 

fi aw{2aw — l) 

kp{x*) 2a-uj—a,{x*) ' 



(37) 



4.2 The Linearized Problem around the Steady State 



In order to investigate local linear stability, we consider in this section the linearization 
around the positive steady state. We first derive a characteristic equation for the eigenvalue 
problem. The signs of the real parts of these eigenvalues give stability (if they all are 
negative) or instability (if there exists one with positive real part). To emphasize our main 
point, which is that stability as well as instability of the positive steady state can take 
place for a suitable choice of model parameters, we shall focus on some simpler cases where 
stability analysis is more transparent. 
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The Characteristic Equation in the General Case 

We denote by {w,u,v) the steady state solution to Equations ( |30|) - (!33|) . Positivity of the 
considered steady state yields a{v) = 0. The linearized problem reads 

= '^{,},.vit), (38) 

dg 

dtu{x,t) + d^[g{x,v)u{x,t)] + d^[—{x,v)u{x)]v(t) = p{x)u{x,t), (39) 

u{0,t) = w{t), (40) 

j^vit)=gix*,v)uix*,t) + ^{x*,v)uix*)v{t) - /it;(t), (41) 

where w, u, and v denote now the deviation of the solution from the steady state. Setting 
w(t) = We^^, u{x,t) = f/(x)e^* and v(t) = Ve^* we obtain the eigenvalue problem of the 
form 

\W = ^(v)wV, (42) 
dv 

XU{x) + d^[g{x,v)U{x)]+d4^ix,v)uix)]V = p{x)U{x), (43) 

f/(0) = W, (44) 

\V = g{x\v)U{x*) + ^(x\v)uix*)V - j^V. (45) 

ov 

Defining an auxiliary function f{x) such that 

f{x)J^'' = -d4^{x,v)u{x)], 



we obtain 



) 



g{x,v)U{x) = giO,v)UiO)e +Ve o^^ / /(s)eo ^'"'^ ds. 



Hence, using ( H2l) and ( l44l) leads to 

/ da w f f -T^d^ \ - f ^f^rfs 



We insert this expression in Equation fl45|) and obtain the characteristic equation 

da f da w ""f I -r^da \ - f ^^^fe^ds 

A + /i-^(x*,^)t.(x*)= U(0,^;)^(t;)^+ //(.)e^^ rfs e . (46) 



dv \ dv \ 
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The Simplest Case: g independent of v 

We first focus on tlie simplest case when the maturation rate g{x,v) = g{x) does not 
depend on f . In other words, the feedback loop only affects stem cells. Although this 
case is very restrictive compared to the original discrete model, since it does not include 
relation (|6]), it is an illustrative example of a possible general behavior. Instability and 
appearance of the oscillations in this model suggest that regulation of the processes solely 
by the stem cell level is not enough to stabilize the system. Moreover, regulatory feedback 
between mature cells and progenitor cells has a stabilizing effect and is essential for efficient 
regulation of the process. 

Since we have / = 0, combining Equation ( 136|) with Equation (146|) . we arrive at 

dcy f 1 

X'^ + fiX = fiv-^{v)e-^^, T= -^ds > 0. (47) 



dv J g{s 



The relationship is identical with the characteristic equation of a delay differential system. 
Indeed, problem f|T8l) - fl2Tl) can be reformulated as a delay differential system. We obtain 
the following result. 

Proposition 4.2 Assume that Equations (124|) - (126|) hold, a(0) > 0, and g is independent 
of V. Consider the steady state {u,v,w) given in Lemma 4-l\ Then, 



(i) fori < TV \-£j{v)\ < ^, the system undergoes a Hopf bifurcation for a single value yUo > 
of the parameter /i. Therefore the steady state can be either locally stable or unstable, 
(a) Further bifurcations also occur for tv 1^(^)1 > 2/c7r + | and k > 1 for at least one 
value /ifc > 0. 

Because in the special case at hand, the system can be reduced to a delay differential 
equation, the linearised stability implies the stability of the nonlinear system, which then 
undergoes a Hopf bifurcation for certain values of the parameters (see [H [16] and the ref- 
erences therein). 

Proof, (i) In order to identify the parameter values for which the bifurcation occurs, we 
look for purely imaginary solutions X = iou with cu G M. We obtain the two following 
relations ^ ^ 

= jjv \ — {v)\ cos(™), TLO = TV \ — {v)\ sin(™). 
dv dv 

By symmetry, we only consider u > 0. The second relation gives a single value tuq G (0, |) 

as soon as 1 < tv \^(v)\ < ^. We can enforce the first relation for a single /i, because 

cos(™o) > 0. This proves statement (i). (ii) For tv \^{v)\ > 2kiT + | and A; > 1, 

the equation tlo = tv \^{v)\sm{Tu) also has a root TUJk G {2k7i,2k'K + |, for which 

cos(rcjfc) > and thus we can find again a fik for which the first equation is satisfied. But 

there might be multiple compatible crossings and several bifurcations are possible. □ 



We now proceed numerically using, for instance Matlab's device DDE BIFTOOL. We 
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check that for the values = r = 1, we get stabihty for jjiv^iy) = —1 and instabihty for 
^iv^iv) = -2. 

This proposition as well as numerical simulations (see Figures H] and |5]) show that instability 
occurs through a Hopf bifurcation, and that regular oscillations appear. 



A Case Motivated by the Discrete Model 

In this section we will study more closely the case given by the relations We can use 
the values of the steady state computed in Equations (I34l) - (l36l) . keeping g{x,v) and p{x) 
fully general. It implies, denoting g{x) = g{x,v) 

da 2ka^p^ kp^ _ g{x*) f P(0 

— — (W) = —-. ; — r- = , W = U[X )exp< — / -^-rdf 

We can show (see the Appendix for detailed calculations) that Equation can be written 

as 

A + ^=-(2a„-l)(^— (-- + l)+y^(..„)-^e. ,(48) 



X 



Case of a{v) = Pw{j^ ~ 1) ^ind g independent of 

Since g{v) is now independent of x, we have that for all x 

g{x, v) = g{0, v)=p^- a{v) = 2p^{l - 



Wkv' 



In particular g[v) = p^ so ^(f) = Pw^- We now substitute g{x,v) = g{v) in Equation 
and obtain 

X* * 

2a^lfp^ fX-p{x)I^ds \^J^ds 

\ + fi = fi — — + 1+ / eo dx]e 

2au, \ A J Pw J 



We compute the first part of the integral term: J — e" ''^ dx = e?™^ — 1, and so 



X* 

\^ ^ ^ T\fP^ > f P^^") — ^ - — ^* UO\ 

A H = (2a^ — 1)1—+ / CP™ dx \ e . (49) 



2au] 2aw V A J p 



Proposition 4.3 Let a{v) he defined by ([6]) with > |, and {u,v,w) be defined by 
Equations fl34l) - fl36|) the unique steady state solution of System pSjl -f l2T|) . If the maturation 
rate g{x,v) is independent of the maturity of the cell x and if the proliferation rate p is 
constant, then the steady state {u, v, w) is locally linearly stable. For a non- decreasing 
proliferation rate, instability may appear. 
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We treat a case of non- decreasing proliferation rate because it is the most biologically 
relevant; however instability may appear for other cases, and is even easier to exhibit, as 
the proof (postponed to the Appendix) shows. Figures [2] and [3] below illustrate a case of 
instability with a nondecreasing proliferation rate. 



5 Numerical Simulations 

In this section we illustrate our theoretical results with a number of numerical simulations. 
We start with a description of our numerical methods. 



5.1 The Numerical Scheme 

We build a simple numerical scheme for System (fT8|) -( l2T]) . We discretize the problem on 
a grid regular in space and adaptive in time. We denote by At^ = t^~^^ — the time step 
between time t^^^ and time t^, by Ax = x*/I the spatial step, where / denotes the number 
of points: Xi = iAx, ^ i ^ I. 

We use an explicit upwind finite volume method for u 



X. 1 



For time discretization, we use a marching technique. At each time t'', we choose the time 
step At^ so as to satisfy the largest possible CFL stability criterion 

At'' 

so that 

A k 

Ar 



Maxxg{x, v'') 

In order to avoid a vanishing time step, it is necessary here to suppose g G L°°. Also, 
more efficient schemes (of WENO type for instance, see p4l [22]) could be used to capture 
discontinuities of g. 

The algorithm is the following: 

• Initialization We use the initial data 



X. 1 
'-2 



• From t'^ to 
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We calculate = a{v'') and define w''^^ = (1 + ^t^a^)w'^. 
We calculate At*^ = ^axtgixr^vk) and define t^+^ = t'' + At^ 
For a boundary condition at z = 0, we define Uq'^^ = w'^^^. 
We define wf"^^ by the following scheme 

-^t^ + A^ = 



We define v'''^^ by 

— — = g[xi, v'juj - f^v ' ^ 

and the term yU?;'^+^in the right hand side is discretized implicitely for stability. 
The reason for the choice of u'} instead of u'}'^^ in the right-hand side of this last 
scheme is due to cell number balance considerations as shown below. 



Cell number balance. From Equation ( fT8l) - (j2Tl) we have obtained the cell number 
balance ( l23l) . We check the equivalent discrete mass balance: 



^7^ + 



At'' ^ Af" At'' 

j=0 



(a^ + g{xo,v''))w'' + ^pjU^Ax - ^iv^^^. 

i=0 



5.2 Numerical Simulations 

First we compare results of the numerical simulations of the discrete and the continuous 
models. To do so, we depart from the discrete values of parameters given in [25]. The 
notations are those of System (I7j)-(l9]), with gi{s, ui) = 2[1 — ai{s)]piUi, pi independent of s, 
di = for i < n and aj(s) = jr^^- It corresponds to the model 1 studied in [T71 [T8] . 



Parameter 


Value 


Parameter 


Value 


Parameter 


Value 




0.77 


Pi 


2.1510-3 day-i 


ds 


0.6925 day-i 


0-2 


0.7689 


P2 


11.2li0-3 day-i 


k 


12.8.10-1° 


0-3 


0.7359 




5.6610-=^ day-i 






04 


0.7678 


Pa 


0.1586 day-i 






05 


0.154 


P5 


0.32 day-i 






Oq 


0.11 


Pe 


0.7 day-i 






a-j 


0.605 


P7 


1 day-i 
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To make comparison easier, for the continuous maturation model, we replace the interval 
[0, X*] by the interval [1, 7] (7 being the number of maturing steps in the discrete model) and 
we define a{x) and p{x) based on parameters in Table [52] by piecewise linear continuous 
functions with values aj and pi at x = i. We take them along a regular grid to obtain 
approximations of a{x) and p{x). 

In Figured], the results of the discrete model are identical with the ones of the continuous 
model if the grid is equal to X = [1,2, ...7] (case / = 6). If the grid becomes finer, we 
observe a slower convergence toward the steady state together with an increase of the 
relative importance of the stem cell population. Though unrealistic from a biological 
viewpoint, it was expected by the derivation of the continuous model from the discrete 
one. It shows that the analogy between the two models is limited. They exhibit different 
quantitative properties (see [18] for a study of the discrete model properties), as well as 
conditions for nontrivial steady state. Moreover, we see that the typical parameter sizes 
have to be adapted. Indeed, the time evolution is much too slow compared to experimental 
data. 

Let us now focus on the stability and instability properties, in order to illustrate the 
theoretical results of Propositions 14.21 and 14. 3[ 

Figures [2] and [3] are an illustration of the instability case stated in Proposition 14. 3[ 
Here, we took a maturity interval [0, X*] with X* = 50, and a proliferation rate p{x) = 
Pu, + Bxx>Y* with pw = 30, Y* = 20, B = 50. We keep a{x) constant equal to = 0.75 
and k = 1.28x10"^ as in the discrete case. The maturation speed g{x) is given by Equation 
([6]). We see that the destabilization is very slow, and our example is very unrealistic, since 
the stem cell population level is tiny. 

In Figures H] and [S] we illustrate instability in the case of Proposition 14.21 We have 
taken here X* = 1, constant proliferation rate p{x) = Pw = Q and maturation rate g{x) = 1 
and a{v) = (^ - l)p^ with = 0.75 and k = 1.28x10"^ 

6 Final Remarks 

In this paper we have developed a structured population model of cell differentiation and 
self-renewal with a nonlinear regulatory feedback between the level of mature cells and the 
rate of the maturation process. We showed that perturbations in the regulatory mechanism 
may lead to the destabilization of the positive steady state, which corresponds to the 
healthy state of the tissue. In particular, we showed that the regulation of stem cells self- 
renewal is not sufficient for stability of the system and the lack of the regulation on the 
level of progenitor cells may lead to the persistent oscillations. This and other stability 
results suggest how imbalanced regulation of cell self-renewal and differentiation may lead 
to the destabilization of the system, which is observed during development of some cancers, 
such as leukemias. The model developed in this paper is rather general and, after adjusting 
it to specific biological assumptions, may serve as a tool to explore the role of different 
regulatory mechanisms in the normal and pathological development. 

Comparing the model to its discrete counterpart we addressed the question of the choice 
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of the right class of models, for example discrete compartments versus continuous matu- 
ration, punctuated by division events. We showed that the models may exhibit different 
dynamics. Interestingly, the structure of steady states varies and the discrete compart- 
mental model admits semi-trivial steady states of the form (0, .., 0, Ui, .., which do not 
exist in the continuous differentiation model. 

To understand the difference between the two models, we derived a limit equation for 
the discrete model assuming that a continuum of differentiation stages can be defined. 
The rationale for such assumption is provided by the fact that differentiation is controlled 
by intracellular biochemical processes, which are indeed continuous in time, at least when 
averaged over a large number of cells. Consequently, for the proper time scaling we have to 
assume that commitment and maturation of cell progenitors do not proceed by the division 
clock (one division = one step in the maturation process) but is a continuous process and 
can take place between the divisions. This observation explains the fundamental difference 
between the two models. The structured population model is indeed not a limit of the 
discrete model with the transitions between compartments correlated to the division of the 
cells. However, the models can exhibit exactly the same dynamics for a suitable choice of 
the maturation rate function g. 




Figure 1: Comparison of numerical simulations using different grids on the interval [1,7], 
from J = 6 (7 points, maturity step dx = 1, discrete model) to / = 100. Left: mature 
cells evolution with time. Right: distribution of cell density along the maturation level, 
at steady state. One can see that the model is extremely sensitive to the number of steps 
(even 7 to 10): small numbers seem to be unstable, whereas for large numbers the numerical 
scheme converges. 
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Figure 2: Example of instability, in illustration of Proposition 14.31 Left: evolution of 
mature cells. Right: evolution of stem cells. 
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Figure 3: Same case as in Figure O Left: final distribution of cells according to their 



maturity level. Right: time evolution of y '^ffi!(i^t)\^dT ' measure the trend to a stable 
maturity level distribution. 
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Figure 4: Example of instability, in illustration of Proposition 14.21 Left: evolution of 
mature cells. Right: evolution of stem cells. 
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Figure 5: Same case as in Figure HI Left: final distribution of cells according to their 

maturity level. Right: time evolution of 
maturity level distribution. 



J \f^u{x,t)\2dx 
f\u{x,t)\'^dx ' 



to measure the trend to a stable 



A Appendix: Proofs of the Results in Section 14.2 



A.l The Characteristic Equation in the Case Derived from the 
Discrete Model 

Using definitions given in ([6]), Equation ( 146|) can be rewritten as 



X a 



A + /i - — (x , v)u(x ) = -- — p.a, — ^ — u(x ) + e-'" a(o / f(s)eo 
av \ law A J 



Calculating the term with / we obtain 



r -^da \ - f -r^i 

i 9(0-, f) ,r _ \ „ ■'. g(s,v) 



ds \e 



(50) 



f ^- do r r)n f ^^^"' dc 

/(x)eo''"'' dx= -d,[^{x,v)u{x)]eo '^'''^ dx 



+ I -^{x,v)u{x)d4e° 



\dx — —[x ,v)u[x je'^ +— (0,?;)m(0) 





X* 



f dg X-pix) f^M"^', dg^ , 

+ / —{x,v)u{x) — - — TT-eo dx-—'- 
J ov g[x,v) ov 





eo dx — Tr{^ ) v)u[x je" + j;—pw'w. 



2a„ 
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Inserting it into Equation ( 150|) leads to 

( k g(x*) , f^*p(Sld£r f 9g , s , .X-p(x) S^^jnr'^'', k .\ - 1 shj<i^ 

V 2a^ A V dv g[x) 2a^ ' ) 



Due to the definition of w given by Equation ( l36i) . the first and the last term of the 
right-hand side can be written together as 



\'^Ojw ^ J civ 9\-^) 



e 



Using Equation ( 13511 we obtain 

u(x) — _, . eo = g(x*)u(x*)e _, ,\ ^ eo ^ ^ . 

g{x) g[xY 

So finally, using Equations (j2S]) and and taking g{x*)u{x*) = ^{2aw — 1), we obtain 
the expression given by Equation (148|1 . 



Proof of Proposition 14.31 

The local linear stability is equivalent to the fact that all eigenvalues A G C, given by 
solutions of Equation fH9|) . have negative real parts. 
First step: A is a solution of the following equation 

X* 

X^ + CX + D = - — (2a^ - 1)A [ 6(x)e^^("-"*)dx (51) 
2a«, J 



with C = > 0, D = Pwl^'^^r^ > 0, and b(x) = > a non-decreasing function. 

Indeed, with this definition of b{x) we rewrite Equation ( H9|) in the form 



A + ^ = -^ 
2a,„ 2a 1 



-^(2a^-l)(^^ + j{l + h{x))e^^dxy~V. 



Integrating by parts we obtain ^ + / e^^^dx = ^ep»^ . 



Second step: the limiting case is for h{x) = 0, i.e., p independent of x. In this case, the 
eigenvalues are given by 

^ -C ± - AD 
A± — 7, • 
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For — AD > these two eigenvalues are negative. If — AD < 0, they are complex 
conjugated with negative real parts. In any case, the steady state is locally linearly stable 
and the first part of the proposition is proved. 

Third step: In the general case, in order to study the sign of the real part of the eigenvalues 
A, we look for values of the parameters such that X = iu with a; e R. It corresponds to a 
Hopf bifurcation and it leads to: 

X* X* 

—uj'^+iCu+D = —-^:j^{2ayj—l)iuj / cos( — (x— x*))(ix+-^(2a^— / b{x)sin{ — {x—x*))dx. 

^Ojw J Pw '^Q'w J Pw 



Taking the imaginary part of this equation yields, since u ^ : 

X* X* 

C = — —{2ayj — 1) / cos(— — x*))dx = — —{2ay, — 1) / b{^* ~ u) cos(— 

2(Jji; J Pw '^O'W J Pw 



Since h is increasing, h{x* — ■) is decreasing. This leads to 

Min{x*,^TT) 

/b{x* — y) cos{—y)dy > 0. 
Pw 



Indeed, either x* < f in which case it is evident because b{x* — y) cos{^y) > for all 
< X < X* , or X* > ^— and we can write 



b{x*-y) cos{— y)dy >b{x*-'^—) / cos{—y)dy-b{x*-^ — 
Pw 2 u J Pw 2 u 



/ cos( — y)dy = 0. 

J Pw 



^pm. 

2 u} 

To end the proof, let us simply exhibit an example where instability can occur: Let x 
denote the Heaviside function. Defining 

b{x) = BXy*<x<x*, 

we compute explicitely 

2 2aw-l 2a^ ~ ^ -o • / ^ / * *nn , 2a^ - 1 . , * *\\ ^\ 

-uj +1- VPwl^^ = -PwfJ'^ iB sm{ — [x -y ))+p^fi— B[cos{ — [x -y 

^Clw ^O'W ^Oiw Pw ^^w Pw 



It provides two relations 



-(2a. - l)B sin{^{x* - y*)), ^ = ^(l + B{l- cos(^(x* - r)))) ■ 
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We can see that there exist many sets of parameters such that both relations are satisfied. 
Given a^, B, p^, and u such that 



< 1, 



PwC^a^ - 1)B 



we can always find x* — y* such that the first relation is satisfied, and then fix n using the 
second relation. n 
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